Spin-orbit coupling and broken spin degeneracy in multilayer graphene 
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Since the lattices of j4_Bj4-stacked graphene multilayers with an even number of layers, as well as 
that of monolayer graphene, satisfy spatial-inversion symmetry, their electronic bands must be spin 
degenerate in the presence of time-inversion symmetry. In intrinsic monolayer and bilayer graphene, 
when symmetry is not broken by external fields, the only spin-orbit coupling present at low energy 
near the corner of the Brillouin zone is the Kane-Mele term, that opens a bulk energy gap but does 
not break the spin degeneracy of the energy bands [C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 
226801 (2005)]. However, spin splitting is allowed in multilayers with an odd number of layers (> 3) 
because their lattices do not satisfy spatial inversion symmetry. We show that, in trilayer graphene, 
in addition to the Kane-Mele term, there is a second type of intrinsic spin-orbit coupling present at 
low energy near the corner of the Brillouin zone. It introduces a Zeeman-like spin splitting of the 
energy bands at each valley, with an opposite sign of the effective magnetic field in the two valleys. 
We estimate the magnitude of the effective field to be ~ 2T. 

PACS numbers: 73.22.Pr 81.05.ue,73.43.Cd. 



Interest in graphene [l|, [1| stems partly from the pres- 
ence of two Dirac cones in the low-energy electronic band 
structure, each of them supporting chiral quasiparticles. 
The origin of the Dirac-like spectrum lies in the fact 
that the honeycomb lattice of graphene has two inequiv- 
alent atomic sites, called A and B, and the amplitude 
of the electronic wave function on them acts as an addi- 
tional degree of freedom known as pseudospin. An 'up' 
component of pseudospin pointing perpendicular to the 
graphene plane would correspond to electronic density 
solely on the A sublattice sites, whereas 'down' pseu- 
dospin corresponds to density on the B sublattice. In 
practice, electronic density is usually shared equally be- 
tween the sublattices so that the pseudospin is a linear 
combination of 'up' and 'down', and it lies in the plane of 
the graphene sheet. There are two Dirac cones, centered 
at inequivalent corners of the Brillouin zone which are 
denoted and and are also referred to as valleys. 
The valleys introduce another spin-like degree of freedom 
into the electronic Hamiltonian. 

The band structure described above is modified when 
the spin of the electron is taken into account. Kane and 
Mele 3 introduced a spin-orbit coupling term that exists 
at the center of each valley, 

Ukm = aU^a^^Sz , (1) 

where a is a parameter and Pauli matrices 11^, cr^ , Sz act 
in K+/K^ valley, A/B sublattice, and f/i spin space, 
respectively. It satisfies time and spatial inversion sym- 
metry requirements by coupling the out-of-plane compo- 
nent of electronic spin with the out-of-plane component 
of pseudospin. This does not break the spin and valley 
degeneracy of the energy bands, but it changes the bal- 
ance between the A and B sublattices, opening a band 
gap and realizing a new, topological state of matter, a 
quantum spin Hall insulator [3(. 

In this paper, we explain how time and spatial- 



inversion symmetry influence spin and valley degeneracy 
in ^i3A-stacked (Bernal) multilayer graphene composed 
of N layers. In the presence of time-reversal symmetry, 
when the atomic lattice satisfies spatial inversion symme- 
try, and such symmetry is not broken by external fields, 
the electronic bands must be spin degenerate. This holds 
for monolayer and bilayer graphene where the only intrin- 
sic spin-orbit coupling term ^-'3] present at low-energy 
near the center of the valley is of the Kane-Mele type, 
Eq. ([T]). However, the lattices of multilayers with odd 
N (TV > 3) do not satisfy spatial inversion symmetry so 
that spin splitting of the energy bands is allowed. We 
show that, in trilayer graphene, in addition to the Kane- 
Mele term, Eq. ([T]), there is a second type of intrinsic 
spin-orbit coupling present at low-energy near the center 
of the valley, 

-Haba = PIlzSz. (2) 

It couples the out-of-plane component of electronic spin 
with the out-of-plane component of the valley 'spin' de- 
gree of freedom. This introduces a Zeeman-like spin split- 
ting of the energy bands at each valley, of magnitude 
proportional to the parameter /3, with an opposite sign 
of the effective Zccman field in the two valleys. 

We begin by briefly describing how time and spa- 
tial symmetries impose constraints on the spectra of 
graphene multilayers, before discussing the particular ex- 
amples of bilayer and trilayer graphene. We consider 
intrinsic spin-orbit terms that exist in samples in the ab- 
sence of external electric or magnetic fields. For example, 
we neglect the influence of a transverse electric field that 
breaks reflection symmetry and produces an additional, 
Rashba-like term 0-11, 0|- Time reversal symmetry re- 
lates energy bands with opposite spin components {f, ],} 
and momenta k and — k: e-|-(k) — e4,(— k). In monolayer 
graphene 3, and in multilayers with even N, includ- 
ing bilayers [l0l - [l3j . the lattice. Fig. [T](a) and (b), obeys 
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FIG. 1: (a) Plan view of the lattice of a single layer of 
graphene showing inequivalent A and B sites, (b) Schematic 
side view of the unit cell of bilayer graphene and (c) of Bernal- 
stacked trilayer graphene, showing the inversion center for 
bilayer and the mirror plane for trilayer. (d) Sketch of the 
trilayer energy bands in the vicinity of the valley center, us- 
ing Eqs. (O and (|B)), where we use a™ = «(, = &m = fet. = 0, 
am = ab = /3m = —f5b = a. with the value a = 71 /20 chosen 
for illustrative purposes. 



spatial inversion symmetry {x,y,z) — > (— x, — y, — z), 
which relates energy bands with the same spin compo- 
nent (as spin is an axial vector) and opposite momenta: 
e-|-(k) — e-|-(— k). Combining time and spatial inver- 
sion symmetry guarantees spin degeneracy of the energy 
bands, e-|'(k) = e4.(k), even in the presence of finite spin- 
orbit coupling. The situation is different in multilayers 
with odd N, starting from trilayer graphene (iV = 3), 
Fig. [IJc), because their lattices do not satisfy spatial in- 
version symmetry (l^ - [l6j . When time reversal symme- 
try holds, the energy of bands with opposite spin and 
momenta must still be equal, e-|-(k) = e4.(— k), but spin 
splitting of the energy bands at the same momentum, 
et(k) ^ e^(k), is allowed, FigM^)- 

Bilayer graphene: In multilayers with an even num- 
ber of layers, the lattice is symmetric with respect to 
spatial inversion symmetry P [(x, y, z) — > {—x, —y, —z)], 
because the point group of the lattice Dsd [iJ, [l3 
({£■, 2C3, 3C2, i, 25*6, 3(Td}) can be regarded as a direct 
product of group D3 ({i?, 2C3, 3C2}) with the inversion 
group Ci {{E,i}). Here, we focus on the particular case 
of bilayer graphene which consists of two cou- 

pled layers of carbon atoms arranged on a honeycomb 
lattice with inequivalent sites {Al, Bl} and {A2, B2} on 
the bottom and top layers, respectively. Fig. (Hb). The 
low-energy tight-binding model takes into account a sin- 
gle Pz orbital per site, with parameter 70 taking into 
account nearest-neighbor hopping, Al-Bl and A2-B2, 
within each layer. The layers are arranged so that sites 



Bl and A2 are directly below and above each other, con- 
nected by interlayer coupling 71 ~ 0.4eV. Then, the band 
structure consists of four bands. Two of them, corre- 
sponding to the coupled orbitals on the Bl and A2 sites, 
are split away from zero energy by ±71. The other two 
bands, formed by effective hopping between Al and B2 
sites, are parabolic and touching near zero energy. 

The parabolic low-energy bands may be described us- 
ing a basis of electronic wave functions with components 
ipi.^.s where i = {Al, B2} labels the relevant lattice sites, 
s — {tii} spin components, and ^ = ±1 distinguishes 
between the valleys located at — ^(|7ra~^,0) with 
lattice constant a. Neglecting spin, the low-energy bands 
are approximately described by the effective Hamiltonian 
[HI = -(«V7i)^o (Hoa, [pl~pl]+2Uzayp,py), 

where we use direct products of Pauli matrices 
o'x,y,z,o'Q = 1 acting in the A1/B2 sublattice space, 
^x,y,z, Ho = 1 acting in the valley space, and Sx,y,z, Sq = 
1 acting in the spin space. Here, v = (•\/3/2) ajo/h is the 
monolayer Fermi velocity and p = {px,Py) is the momen- 
tum measured with respect to the center of the valley, 
k = KK^ + p. 

We begin by exploring which spin-orbit terms, if any, 
satisfy time and space symmetries of the lattice. Apart 
from symmetry, we assume that low-energy states are 
composed from orbitals on the Al, B2 sites, so that the 
following argument does not rely on details of a particular 
model. In the above basis, the time reversal operator is 
T ~ \lx<^QSyJC where JC represents complex conjugation, 
and translation by a distance equal to the lattice constant 
a along the x-axis is given by e^^^^^/^ . Generators of the 
point group of the lattice Dy,^ are an active rotation anti- 
clockwise by angle 27r/3, C3 = _e^'^n.<T,/3g-i7rS,/3^ j.^^^^, 

tion by angle tt about the x axis, C2 = —iHoaxSx, and 
spatial inversion, i = Ilx<yxSo. We focus on the center of 
the valley, p = |p| = 0, and consider all possible terms 
describing coupling between valley, sublattice and spin of 
the form HiCTjSk where i — {{),x,y,z}, j = {0,x,y, z}, 
k — {0, X, y, z}. Then, we determine which of them is in- 
variant with respect to the time and space symmetries of 
the bilayer lattice. Except for a constant term HoctqS'o, 
there is only one, which is a term of the Kane-Mele type, 
Eq. ID), [^-Q. As in a monolayer ^, it leads to the open- 
ing of a gap at the center of the valley, but doesn't break 
spin and valley degeneracy. 

Trilayer graphene: In odd-A'^ multilayers with A^ > 3, 
the point group {{E, 2C3, 3C^, a?,, 2^3, Sa,}) 
can be regarded as a direct product of the group D3 
({£:, 2C3,3C^}) with the reflection group {{E,ah}). 
With respect to the even-A^ multilayers, spatial inversion 
is replaced by mirror reflection ah [{x, y, z) — {x, y, —z)]. 
Trilayer graphene [H, [ij, [T6l - [23 | consists of three cou- 
pled layers of carbon atoms arranged on a honeycomb 
lattice with inequivalent sites {A1,B1}, {^2,B2} and 
{A3, Bi} on the bottom, middle and top layers, respec- 
tively, Fig.[ljc). The simple tight-binding model using a 
single Pz orbital per atom includes coupling 70 between 
nearest-neighbor sites within each layer, Al-Bl, A2-B2 



3 



and A3-B3, as well as interlayer coupling 71 between sites 
Bl and A2, and A2 and B3, that are directly below and 
above each other. 

Linear combinations of atomic orbitals may be sep- 
arated [1^ [1^ into a pair that are odd with respect 
to mirror reflection symmetry cr/,,, (f>^"^\ = [V'Ai.j.s ~ 



Hamiltonian. We consider the two low-energy bands 
that touch near zero energy, related to orbitals cfif''^ ^ , 
^■il s- ^'^^ ^^^^ basis, T = IlxCr^^ SylC, translation is 

and ail = ^'i^o^'o''^ Sz- The invariant terms at the center 



■0^3,4, s]/V2, = [V'Bi.f.s - V'B3,e,s]/V2, and four of the valley are 



that are even with respect to ah, 4'i'\s ~ [V'Ai.^.s 



(b) 



3,5, s 



bpBi,^,s+ipB3,^,s]/V^- Written in a basis of states 0-™{''\ 
the Hamiltonian of trilayer graphene separates into two 
parts The first has a basis of odd orbitals ^i™^^, 

02™ ''s ' ^^'^ ^ form analogous to that of monolayer 

graphene. 



(3) 



where matrices crQ™'^^^ act in the 0^™^ , 02™'' sublattice 
space. The second part of the Hamiltonian has a basis 



of four even orbitals 



,(b) 



^ J, , resulting in four parabolic 
bilayerlike bands. As in a bilayer, two bands touch near 
zero energy. They are related to orbitals ^^''^ ^ , (/ij''^ ^ and 
may be described by an effective quadratic Hamiltonian 
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(0) 
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[Hoai^) (pI -pI) + 2IVzaf^p,py\ , (4) 



where matrices a, 



(b) 



act in the 



sublattice 



space. 

An analysis of the transformation properties of the tri- 
layer Hamiltonian under symmetries of the lattice is rela- 
tively simple in the basis of orbitals 0-™{*^ because time 
and spatial symmetry operators do not mix the mono- 
layerlike and bilayerlike parts of the Hamiltonian. We 
therefore consider them separately. In the monolayerlike 
part, the time reversal operator is T = Tl^aj^^ SylC, and 
translation by a distance equal to the lattice constant a 
along the x-axis is given by e*'*^'^^/'^. Generators of the 
group 1)3/1 are an active rotation anti-clockwise by angle 

2^3, C3 = e-^'^^^/3[(l-ai™))/2+e-*2.n^/3(l+ai'"V2], 
rotation by angle tt about the y axis, C'2 = iHxaj^"^ Sy, 

and mirror reflection symmetry, a^, — iHoa^'^ Sz- 

We focus on the center of the valley, p — 0, and 
consider all possible terms describing coupling between 
valley, sublattice and spin of the form IliajSk where 
i = {0,x,y,z}, j = {0,x,y,z}, k = {0,x,y,z}. Then, 
we determine which of them is invariant with respect to 
the time and space symmetries of the lattice: 

ni\^ = a™Ho4"')5o + 6™Hoai'")S'o 

+ a„H,a(™)5, -I- /3„H,a^"^5, , (5) 

The first line of Eq. ([S]) contains spin-independent terms, 
the second line contains spin-orbit terms. We repeat 
the symmetry analysis for the bilayerlike part of the 
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(1) 



Hoa^'^S-o + bbHoai^^^o 
Uzai'^Sz+mz^'^Sz, 



OLb 



(6) 



and they are analogous to the terms in the monolayerlike 
part Eq. 

The monolayerlike ["Hm' +'Hni] and bilayerlike [h^^^ + 
'H\^^] parts of the Hamiltonian produce superimposed 
monolayerlike tm and bilayerlike bands, 

Em = a,n + siPra ± \/ (6m + S^Om)^ + W^P^ , (7) 

£6 = a6 + se/3f, ±^{h + s^abf + i^V/(27?) , (8) 

where s = ±1 denotes different spin components, ^ = 
±1 different valleys, and we consider {|em/6|,wp} ^ 7i. 
Fig. [lid) shows a plot of the trilayer energy bands, using 
Eqs. ([7]) and ([5]), where we use a™ = ai, = b„i = bh = 0, 
ctm — OLb — Pm = —Pb = OL and the value a = 71 /20. 

The terms proportional to the parameters Om and ab 
in Eqs. ^ and ^ produce constant energy shifts of the 
monolayerlike and bilayerlike spectra with respect to each 
other. The term proportional to hm {bb) is absent in 
a monolayer (bilayer), but allowed in trilayers owing to 
the absence of spatial inversion symmetry. It breaks the 

01™'' 702™'' (0i''V02'''') sublattice symmetry and opens a 
gap between the monolayer (bilayer) bands, but doesn't 
generally break valley symmetry The Kane-Mele 

spin-orbit term, responsible for factors am/6, has been 
discussed in detail elsewhere [^-Q- Here, we note that, 
in trilayers, it can conspire with the 6m/6 parameters to 
break spin degeneracy even at zero magnetic field. Pa- 
rameters Pm/b arise from the spin-orbit term unique to 
odd-A'^ multilayers and they describe a Zeeman-like spin 
splitting of the energy bands at each valley, with an op- 
posite sign of the effective Zeeman field in the two valleys. 

The Hamiltonians Eqs. ^ and ([6]) were obtained 
using symmetry arguments, but it is possible to relate 
the phenomenological parameters to tight-binding pa- 
rameters. For example, for the spin-independent part 
of the Hamiltonian, parameters am, 6m, ifc, may be 
related to next-nearest-layer coupling parameters 72 
tween Al and A3) and 75 (between Bl and B3) as 
a„i = -(72 + 75)/2, bm = -(72 - 75)72, ftb = 7272, 
bb = 7272. 

While it is generally accepted that the magnitude a 
of the intrinsic spin-orbit couplixig near the center of the 
valleys in graphene is small [1, |23[ . its value has been 
the subject of theoretical debateji, i, 0, i, [lUli]. It 
was recently proposed j26|] that coupling between the Pz 
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orbitals and d orbitals gives the dominant contribution 
to spin-orbit coupling in monolayers, because Pz orbitals 
are not orthogonal to d^z and dyz orbitals, yielding a ~ 
O.OlmeV. Very recently, it has been suggested 0, !3| that 
spin-orbit coupling in bilayer graphene can be relatively 
large due to coupling between pz orbitals and a bands 
on different layers (that are not orthogonal), giving a 
contribution to the spin-orbit coupling a ~ O.lmeV. 

Using the tight-binding model, we generalize the esti- 
mate of the magnitude of spin-orbit coupling in mono- 
layers d, d, [13) and bilayers 0, Q to demonstrate that 
both types of intrinsic spin-orbit coupling are likely to 
be as large in trilayers as the Kane-Mele term in bilay- 
ers: |/3m/f)| ^ \oim/b\ ^ O.lmeV. We take into account 
coupling of the pz and a orbitals, by writing the tight- 
binding Hamiltonian including the Hamiltonian of the Pz 
orbitals and of the cr orbitals H„ as 



H 



(9) 



where matrix V describes coupling between the pz and 
a orbitals. Here, is a 12 by 12 matrix (describing a 
spin up and spin down pz orbital per site, and six sites 
in the trilayer unit cell) and H„ is a 36 by 36 matrix (six 
orbitals per site taking into account two spins and s, p,^. 
and py , and six sites in the trilayer unit cell) . 

Block H„ contains inter-atomic matrix elements that 
are written using Slater and Koster matrix elements [27| . 
For example, at the valley center, the matrix element 
between a Px orbital on an Aj site and a py orbital on 
an adjacent Bj site within the same layer {j — 1,2,3) 
is ((p^^\H\p^^) = --iiSypp/A. Matrix V in Eq. © de- 



scribes the coupling between pz orbitals and a orbitals. 
We take into account skew interlayer coupling between 
orbitals on sites Al^ B2 and A3, that are not coupled 
by 7i (such coupling, if it is only between pz orbitals, is 
usually denoted 73). For example, the matrix elements 
between a pz orbital on the Al site and a px or py or- 
bital on the B2 site are (pf^\H\pl^'^) = 3i£,Vpp/2 and 
{pf^\H\p^^) = -Wpp/2 exactly at the valley center. For 
an A3 site instead of ^41, these matrix elements acquire 
an additional minus sign. 

For the spin-orbit interaction, we include intra-atomic 
coupling H = AqL • S between the p orbitals. This in- 
troduces spin-dependent matrix elements within for 
the same atomic site j such as {p^. ^\H\py ^) — — zAq, 

(Px l\-^\Py i) ^ *^o- "^^"^ effective Hamiltonian de- 
scribing the Pz orbitals at zero energy may be writ- 
ten using second-order perturbation theory as H^s ~ 
Htt ~ VH^^V^. Doing so, we recover the intrinsic spin- 
orbit terms, Eqs. (O and ([S]), where am — ctb = Pm — 
~l3b — 2Ao(V^p/lpp)^. The magnitude of this estimate 
is the same as that for the Kane-Mele term in bilayer 
graphene 0. Using Vp^/Ypp ^ 0.1 and Aq ^ lOmeV 
yields lam/hl = |/3m/fe| ~ O.lmeV. We compare the spin 
sphtting of the bands, 2|/3„j/fc|, to the Zeeman energy 
gUsB in a real external field B. Using g = 2 gives an 
effective field of i? ^ 2T, which is comparable to the 
real fields required to observe Zeeman-split conductance 
fluctuations in graphene [23|. 
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